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I. Introduction 

(U) In combustion computations, it is well known that the predictions 
of chemical reaction rates (the source terms in the species conservation equa- 
tions) are poor if conventional turbulence models are used. The main dif- 
ficulty lies in the fact that the reaction rate is highly non-linear, and the 
use of averaged temperature, pressure and density produces excessively large 
errors. Moment closure models for the source terms have attained only lim- 
ited success. The probability density function (pdf) method seems to be 
t e only alternative at the present time that uses local instantaneous val-. 
ues of the temperature, density, etc., in predicting chemical reaction rates 

and thus is the only viable approach for more accurate turbulent combustion 
calculations. 

, Th ** act tllat tlle I 3cit equation has a very large dimensionality ren- 

ers finite difference schemes extremely demanding on computer memories 
and thus impractical, if not entirely impossible. A logical alternative is the 
Monte Carlo scheme, which has been used extensively in statistical physics. 
The evolution equations for the joint pdf of the velocity and species mass 
fraction have been successfully solved using Monte Carlo schemes see e - 
Pope[l]. However, since CFD has reached a certain degree of maturity °aa’ 
well as acceptance, it seems, at least from the stand point of practical appli- 
cations. that the use of a combined CFD and Monte Carlo scheme is more 
beneficial. Therefore, in the present study a scheme is chosen that uses a 
conventional CFD flow solver in calculating the flowfield properties such as 
velocity pressure, etc., while the chemical reaction part is solved using a 
Monte Carlo scheme. 

(U) A combined CFD- Monte Carlo computer algorithm has been de- 
veloped recently. As a first calibration of the Monte Carlo solver recently 
developed in this work, the discharge of a heated turbulent plane jet into 
quiescent air is studied. Experimental data for this problem shows that when 
the temperature difference between the jet and the surrounding air is small, 
buoyancy effects can be neglected and the temperature can be treated as a 
passive scalar. The fact that jet flows have a self-similar solution lends con- 
venience m the modeling study. Furthermore, the existence of experimental 
ata for turbulent shear stress and temperature variance (temperature fluc- 
tuation) make the case ideal for the testing of pdf models wherein these 
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values can be directly evaluated. 

(U) The following section presents the methodologies used in this study, 
and a discussion on the numerical results and their comparison with experi- 
mental data follows. 


II. Method 

2.1 Governing equations 

(U) The Favre averaged momentum and species transport equations can 
be written as 

pd t Ui + pujdjUi = - dip + vdjTij — dju"u" 

pd t Y\ + p&jdjYi = pdjiDdjYi - V'"u") + d^Ddjpd’Yn + pw t 

where u, is the velocity and Y t is the mass fraction. The major problem en- 
countered in the numerical study of turbulent combustion lies in the modeling 
of the chemical source term, pw^ known as the chemistry closure problem. 
The Source term w t is an exponential function of temperature, and it is well 
know that the use of averaged temperature, T 1 , in evaluating pw{ can cause 
egregious errors. Therefore the accurate prediction of turbulent combustion 
using conventional turbulence models becomes very difficult. This motivated 
the use of pdf methods. If the pdf, P. is given, then the mean of the source 
term can be evaluated exactly: 

pwi = J-.JpwiM Y n ,T,p)P(Y u ...,Y n ,T,p)dY 1 ...dY n dTdp 

To solve for the pdf, we need the following: 

2.2 Evolution equation for the pdf 

(U) The evolution equation for the probability density function of the 
mass fractions, temperature, and density, P(Y\ y p), can be written 
as[i] 

N 

pd t P + pv 0 d a P + p'^2^,{w i (xl; 1 ,...,ip N )P} 

1=1 

= -d a (p < t CW’ t > P) -PEEdU« tiiWk > p) 

.=1 J=1 
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where the terms represent mean convection, chemical reactions, turbulent 
convection, and molecular mixing, respectively; P is the density-weighted 
joint pdf: 

p = pp, 

e is the scalar dissipation: 

= ~pD tjdjY,. 

ipi s are variables such as the mass fractions, density, and temperature, and 
<, > denotes the mathematical expectation ol a function. 

(L) The left hand side of the above equation can be evaluated exactly 
and requires no modeling; the right hand side terms contain the conditional 
average of the Favre velocity fluctuation and the conditional average of the 
scalar dissipation and require modeling. 

2.3 Finite difference solution for the NS equations 

(U) In oidei to simplify the problem so as to concentrate on the study of 
the pdf models, we have confined our numerical procedures to parabolic flows. 
A parabolic NS solver with a k-e turbulence model has been developed. 
The k - e model has been tested by solving standard cases such as flat plate 
boundary layers and free shear layers with satisfactory numerical results. 

2.4 Monte Carlo scheme for the pdf equation 

(U) A grid dependent Monte Carlo scheme has been employed, primarily 
following Chen and Kollmann[5], The task here is to construct an ensemble 
of sample points, each sample has its own distinct properties such as tem- 
perature. mass fraction, etc.: these properties change with time or location 
such that the probability function of the ensemble evolves according to the 
evolution equation. Consequently, the pdf of the ensemble is the desired 
approximate solution for the pdf equation. 

(U) Suppose N samples are assigned to a grid cell in the flow domain, 
the pdf for the ensemble of N samples can then be written as 

P’tt’) = T dj) 

i=i 

where o } is the scalar function value carried by the j th sample, e.g., the mass 
fractions, etc. c 
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(U) In order to have the pdf P m evolve according to the pdf evolution 
equation, we discretize the equation on a given grid and write, for parabolic 
flows: 

P*z+dz,j = &jP*x.j + l "T djP x,j "l" 7 jP xj - 1 

and require 

otj + 3j + 7i = 1 

2.5 Recontamination in a Monte Carlo scheme 

(U) Using a very simple test case of a convection/ diffusion process with 
two scalars, it was found that the previous scheme[5] does not conserve mass 
fractions due to re-contamination. It is found that in order to conserve the 
mass fractions absolutely, one needs to add further restriction to the scheme, 
namely 

<Xj + Ij - + T)+1 (*) 

A new computer algorithm was devised and tested. This algorithm uses a 
few extra arrays to store informations form the previous time level and thus 
eliminates the repetition in the sampling process. With the simple test case 
of two scalars with assumed constant coefficients in the pdf equation, the 
new algorithm is shown to conserve the mass fractions perfectly. Deficiencies 
such as directional bias and re-contamination that were found in the previous 
algorithm are completely eliminated. 

(U) It is not yet known whether one can indeed devise a scheme that 
satisfies relation (*) in general flows, where the coefficients of the pdf equation 
are variable and are calculated from the flow velocities, turbulence time scale, 
etc. 
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ILL Results 

(U) The sketch of a heated plane jet is given on this page. The flow 
domain is divided into 30 cells in the cross direction, and 100 samples are 
assigned to each cell. The numerical results from the present study are 
compared with experimental data in the following figures. 



Figure 1. Sketch of a heated turbulent free jet. 
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(U) The numerical solution for the ensemble averaged pdf’s of the tem- 
perature: y/yo .5 = 0. where y 0 .s is the location where T — 0.5T o . In order to 
interpret the data, the ensemble averaged pdfs are integrated to obtain the 

approximate continuous pdf: P = P*(AT)/AT. The result of this integra- 
tion is shown as vertical bars. The normal distribution is plotted as a dotted 
line for comparison. 


Fig. 2 probability density distribution 



(T- < T> )/ a 
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P(kT)/M 


(U) The numerical solution for the ensemble averaged pdf’s of the tem- 
perature: y/yo.s = 0.44, where y 0 .s is the location where T = 0.5 T 0 . In order 
to interpret the data, the ensemble averaged pdfs are integrated to obtain 
the approximate continuous pdf: P = F*(AT)/AT. The result of this in- 
tegration is shown as vertical bars. The normal distribution is plotted as a 
dotted line for comparison. 


Fig. 3 probability density distribution 
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(U) The computed velocity profile is compared with experimental data 
from Foerthmann[7], The good agreement between the numerical solution 
and experiments indicates that the parabolic flow solver employed in the 
present study is fairly accurate and can serve its purpose. 


m 


Fig. 4 Velocity distribution 




(U) As a calibration for the k - t turbulence model, the turbulent 
shear stress from the present study is compared with measurements given 
by Brandbury[S], and by Gutmark and Wygnanski[9], Except for a slight 
over prediction towards the edge of the jet, the numerical solution generally 
agrees well with the experimental data. This comparison is very crucial for 
pdf calculations because it is an indicator of the accuracy of the numerically 
predicted turbulent kinetic energy, k. and dissipation, e, which are used to 
determine turbulent time scale, t, from t = k/e. An accurate prediction for 
k and c is essential for the correct modeling of molecular mixing in the pdf 
calculation. 


Fig. 5 Turbulent shear stress 
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(U) The standard variance of the temperature, or the temperature fluctu- 
ation, compared with experimental data from Ref. 12 and 13. Because of the 
limited number of sample points (N=100) used in the computation, the nu- 
merical results are scattered; however, these results agree very well with the 
experimental observations. Particularly important is the agreement between 
the predicted and the measured temperature fluctuations, for in case of finite 
rate chemistry calculations, the instantaneous temperature will be used in 
determining the reaction rate, and the correct prediction of the temperature 
distribution will ensure the accuracy in reaction rate calculations. 



y/y.5T 0 
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